Trait Ordination Plant Community Analysis

Trait Ordination Plant Community Analysis

1. traits

   1. above- and belowground

      1. height & SLA
         

   2. acquisitive vs. conservative resource-use strategies

      1. rooting depth

   3. collaborative vs. individual resource-use strategies

      1. fine root proportion

      2. thickness of roots

Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.

Author: Caryn D. Iwanaga

Updated: 02/06/2025

Libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)

Read in Data

cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
  mutate(
    nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
    ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))


# Define a reusable color scale
# custom_colors <- scale_fill_manual(
#   values = c(
#     "control" = rgb(195, 197, 193, maxColorValue = 250),
#     "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
#     "+compost" = rgb(167, 156, 109, maxColorValue = 225)
#   )
# )

custom_colors <- scale_fill_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

# Set a global theme
theme_set(theme_bw())

Make site x species matrix

1. Calculate relative abundances of species for each site
   1. confused how to aggregate abundances because even grouping by plotid (nut_trt and ppt_trt), year, species, and grazing history (block 1&2 – high, block 3&4 – low) and averaged, the sd are as high as 45.96!
      1. proceeded with these high sd to see what happens (ask Lauren)
   2. try how it is when you split everything by each block (more representative?)

site.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  mutate(
    grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
  ) %>% # separate by block (grazing intensities)
  group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

matrix0 <- spread(
  site.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix <- matrix0[, c(5:79)]

# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)

Split matrix by each block as well

site.block.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

site.block.sd.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    pct_cover,
    abundance = mean(pct_cover),
    sd = sd(pct_cover)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.

matrix0.block <- spread(
  site.block.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0.block) <- paste(matrix0.block$nut_trt, matrix0.block$ppt_trt, matrix0.block$yr, matrix0.block$block, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix.block <- matrix0.block[, -c(0:4)]

# relativize data with Bray-Curtis dissimilarity
matrix.block.bray <- decostand(matrix.block, "total")
rownames(matrix.block.bray) <- rownames(matrix0.block)

Run NMDS

Rules of thumb:

• < 0.05 is excellent

• 0.05-0.1 is good

• 0.1-0.2 is fair

• 0.2-0.3 is cause for concern…

• > 0.3 is poor, random

All years – low & high grazing

Stress: fair

0.1966553 

Grouping by year and grazing history had most significant clustering shown.

metaMDS()

nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)

Run 0 stress 0.2020736 
Run 1 stress 0.2037794 
Run 2 stress 0.1973259 
... New best solution
... Procrustes: rmse 0.07504543  max resid 0.2323388 
Run 3 stress 0.1970656 
... New best solution
... Procrustes: rmse 0.02167554  max resid 0.1184388 
Run 4 stress 0.2112285 
Run 5 stress 0.2016301 
Run 6 stress 0.2034469 
Run 7 stress 0.2100941 
Run 8 stress 0.2020249 
Run 9 stress 0.1990567 
Run 10 stress 0.205915 
Run 11 stress 0.2072425 
Run 12 stress 0.2077976 
Run 13 stress 0.196667 
... New best solution
... Procrustes: rmse 0.01874332  max resid 0.1106696 
Run 14 stress 0.2021165 
Run 15 stress 0.2061482 
Run 16 stress 0.1966667 
... New best solution
... Procrustes: rmse 0.0005495779  max resid 0.002877469 
... Similar to previous best
Run 17 stress 0.2158361 
Run 18 stress 0.21463 
Run 19 stress 0.2120927 
Run 20 stress 0.2037336 
*** Best solution repeated 1 times

nmds.comp # pretty not great stress

Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.1966667 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray’ 

stressplot(nmds.comp)

plot(nmds.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))

nmds_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years – blocks

Stress: cause for concern

0.2329306

metaMDS()

nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)

Run 0 stress 0.2329306 
Run 1 stress 0.2385405 
Run 2 stress 0.2608291 
Run 3 stress 0.2618978 
Run 4 stress 0.2625943 
Run 5 stress 0.2369671 
Run 6 stress 0.2329931 
... Procrustes: rmse 0.002596777  max resid 0.01287822 
Run 7 stress 0.2447034 
Run 8 stress 0.2406908 
Run 9 stress 0.2476092 
Run 10 stress 0.2412944 
Run 11 stress 0.232978 
... Procrustes: rmse 0.002270404  max resid 0.012997 
Run 12 stress 0.23922 
Run 13 stress 0.2497086 
Run 14 stress 0.2478848 
Run 15 stress 0.2365873 
Run 16 stress 0.2329684 
... Procrustes: rmse 0.001740474  max resid 0.01298233 
Run 17 stress 0.2520339 
Run 18 stress 0.248187 
Run 19 stress 0.2380139 
Run 20 stress 0.2333533 
... Procrustes: rmse 0.0115742  max resid 0.112582 
Run 21 stress 0.2489426 
Run 22 stress 0.2400835 
Run 23 stress 0.2421075 
Run 24 stress 0.2510936 
Run 25 stress 0.232944 
... Procrustes: rmse 0.001509127  max resid 0.01190131 
*** Best solution was not repeated -- monoMDS stopping criteria:
    24: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

nmds.block.comp

Call:
metaMDS(comm = matrix.block.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.block.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.2329306 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.block.bray’ 

stressplot(nmds.block.comp)

plot(nmds.block.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))

nmds_dat.block <- cbind(nmds1, nmds2) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year - Blocks")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = as.factor(block), color = as.factor(block))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Block")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  color = as.factor(block), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(block)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Block & Year - Blocks") +
  scale_color_manual(
  values = c(
    # "low" = "cyan",
    # "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years - LOW GRAZING ONLY

Stress: fair

0.1869068

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
  filter(grazing_hist == "low")

row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")

matrix.low <- matrix.low[, c(5:79)]

rownames(matrix.low) <- row.low

Warning: Setting row names on a tibble is deprecated.

matrix.bray.low <- decostand(matrix.low, "total")

metaMDS()

nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)

Run 0 stress 0.2000574 
Run 1 stress 0.1960027 
... New best solution
... Procrustes: rmse 0.03759248  max resid 0.136796 
Run 2 stress 0.2299212 
Run 3 stress 0.1908791 
... New best solution
... Procrustes: rmse 0.1110161  max resid 0.2532458 
Run 4 stress 0.1869068 
... New best solution
... Procrustes: rmse 0.1544638  max resid 0.4035849 
Run 5 stress 0.2066946 
Run 6 stress 0.1878648 
Run 7 stress 0.2004585 
Run 8 stress 0.1876855 
Run 9 stress 0.2113763 
Run 10 stress 0.2243075 
Run 11 stress 0.1964757 
Run 12 stress 0.1893384 
Run 13 stress 0.2100567 
Run 14 stress 0.2455646 
Run 15 stress 0.1930573 
Run 16 stress 0.2043372 
Run 17 stress 0.245567 
Run 18 stress 0.2240266 
Run 19 stress 0.1942011 
Run 20 stress 0.1893384 
Run 21 stress 0.2104986 
Run 22 stress 0.1960027 
Run 23 stress 0.22813 
Run 24 stress 0.2593359 
Run 25 stress 0.1908791 
*** Best solution was not repeated -- monoMDS stopping criteria:
    23: stress ratio > sratmax
     2: scale factor of the gradient < sfgrmin

nmds.comp.low # pretty not great stress

Call:
metaMDS(comm = matrix.bray.low, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.low 
Distance: bray 

Dimensions: 2 
Stress:     0.1869068 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 4 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.low’ 

stressplot(nmds.comp.low)

plot(nmds.comp.low, type="t")

nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
  )
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))

nmds_low_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  color = as.factor(yr), shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

# add year ellipses ---- 
ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
    stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

Make each individual year NMDS to see if there’s better stress

2019

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
  filter(yr == 2019)

row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")

matrix.2019 <- matrix.2019[, c(5:79)]

rownames(matrix.2019) <- row.2019

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2019 <- decostand(matrix.2019, "total")

nmds 2019

Stress: Excellent

0.0560178

metaMDS()

nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)

Run 0 stress 0.0560178 
Run 1 stress 0.06223133 
Run 2 stress 0.05984339 
Run 3 stress 0.05905532 
Run 4 stress 0.06384512 
Run 5 stress 0.05851802 
Run 6 stress 0.05852022 
Run 7 stress 0.06942148 
Run 8 stress 0.05654879 
Run 9 stress 0.06942148 
Run 10 stress 0.05654879 
Run 11 stress 0.05877078 
Run 12 stress 0.0662005 
Run 13 stress 0.05982901 
Run 14 stress 0.05957995 
Run 15 stress 0.07018407 
Run 16 stress 0.05652852 
Run 17 stress 0.06074817 
Run 18 stress 0.05652852 
Run 19 stress 0.05654879 
Run 20 stress 0.07335647 
Run 21 stress 0.05851807 
Run 22 stress 0.05652852 
Run 23 stress 0.0735608 
Run 24 stress 0.0761203 
Run 25 stress 0.06228134 
*** Best solution was not repeated -- monoMDS stopping criteria:
    13: stress ratio > sratmax
    12: scale factor of the gradient < sfgrmin

nmds.2019 # pretty not great stress

Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2019 
Distance: bray 

Dimensions: 2 
Stress:     0.0560178 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2019’ 

stressplot(nmds.2019)

plot(nmds.2019, type="t")

nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
  )
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))

nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019 - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019 - Grouping by Precipitation Treatment") 

2020

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
  filter(yr == 2020)

row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")

matrix.2020 <- matrix.2020[, c(5:79)]

rownames(matrix.2020) <- row.2020

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2020 <- decostand(matrix.2020, "total")

nmds 2020

Stress: fair

0.1522495

metaMDS()

nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)

Run 0 stress 0.1524241 
Run 1 stress 0.1584809 
Run 2 stress 0.1584809 
Run 3 stress 0.1688381 
Run 4 stress 0.1625903 
Run 5 stress 0.153917 
Run 6 stress 0.3515938 
Run 7 stress 0.1692911 
Run 8 stress 0.1692911 
Run 9 stress 0.1584809 
Run 10 stress 0.1522495 
... New best solution
... Procrustes: rmse 0.02107749  max resid 0.06346372 
Run 11 stress 0.1584809 
Run 12 stress 0.1522045 
... New best solution
... Procrustes: rmse 0.01183105  max resid 0.0373035 
Run 13 stress 0.1744238 
Run 14 stress 0.1524241 
... Procrustes: rmse 0.02481397  max resid 0.06599247 
Run 15 stress 0.1742461 
Run 16 stress 0.1595062 
Run 17 stress 0.1524241 
... Procrustes: rmse 0.02481007  max resid 0.06599555 
Run 18 stress 0.1625903 
Run 19 stress 0.1522495 
... Procrustes: rmse 0.01183378  max resid 0.03737577 
Run 20 stress 0.1522495 
... Procrustes: rmse 0.0118299  max resid 0.03737025 
Run 21 stress 0.1747777 
Run 22 stress 0.184858 
Run 23 stress 0.1524241 
... Procrustes: rmse 0.02481216  max resid 0.06599381 
Run 24 stress 0.1866314 
Run 25 stress 0.2912948 
*** Best solution was not repeated -- monoMDS stopping criteria:
    20: stress ratio > sratmax
     5: scale factor of the gradient < sfgrmin

nmds.2020

Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2020 
Distance: bray 

Dimensions: 2 
Stress:     0.1522045 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 12 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2020’ 

stressplot(nmds.2020)

plot(nmds.2020, type="t")

nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
  )
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))

nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 1) + 
  labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment") 

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Nutrient Treatment")

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Precipitation Treatment")

2021

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
  filter(yr == 2021)

row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")

matrix.2021 <- matrix.2021[, c(5:79)]

rownames(matrix.2021) <- row.2021

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021 <- decostand(matrix.2021, "total")

2021 - LOW

Stress

0.08067789

# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "low") 

row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")

matrix.2021.1 <- matrix.2021.1[, c(5:79)]

rownames(matrix.2021.1) <- row.2021.1

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")

nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)

Run 0 stress 0.08067789 
Run 1 stress 0.1407365 
Run 2 stress 0.131677 
Run 3 stress 0.08067789 
... New best solution
... Procrustes: rmse 4.448434e-06  max resid 7.224535e-06 
... Similar to previous best
Run 4 stress 0.1393679 
Run 5 stress 0.08067789 
... New best solution
... Procrustes: rmse 1.053765e-06  max resid 1.585547e-06 
... Similar to previous best
Run 6 stress 0.1393365 
Run 7 stress 0.1393365 
Run 8 stress 0.08067789 
... Procrustes: rmse 5.054553e-06  max resid 8.236855e-06 
... Similar to previous best
Run 9 stress 0.08067789 
... Procrustes: rmse 5.780641e-06  max resid 8.023545e-06 
... Similar to previous best
Run 10 stress 0.08067789 
... Procrustes: rmse 4.42502e-06  max resid 6.754172e-06 
... Similar to previous best
Run 11 stress 0.2039408 
Run 12 stress 0.08067789 
... Procrustes: rmse 1.729526e-06  max resid 2.63632e-06 
... Similar to previous best
Run 13 stress 0.08067789 
... New best solution
... Procrustes: rmse 1.098472e-06  max resid 1.772254e-06 
... Similar to previous best
Run 14 stress 0.1316771 
Run 15 stress 0.2639821 
Run 16 stress 0.203941 
Run 17 stress 0.08067789 
... Procrustes: rmse 3.956956e-06  max resid 5.336733e-06 
... Similar to previous best
Run 18 stress 0.08067789 
... Procrustes: rmse 1.630904e-06  max resid 3.249717e-06 
... Similar to previous best
Run 19 stress 0.1203381 
Run 20 stress 0.2039414 
*** Best solution repeated 3 times

nmds.2021.1

Call:
metaMDS(comm = matrix.bray.2021.1, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021.1 
Distance: bray 

Dimensions: 2 
Stress:     0.08067789 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021.1’ 

stressplot(nmds.2021.1)

plot(nmds.2021.1, type="t")

dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))

nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))

nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
#   geom_point(size = 4, stroke = 0.5) +
#   labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")

2021 - HIGH

Stress

0.08065034

# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "high") 

row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")

matrix.2021.2 <- matrix.2021.2[, c(5:79)]

rownames(matrix.2021.2) <- row.2021.2

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")

nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)

Run 0 stress 0.08065034 
Run 1 stress 0.1378308 
Run 2 stress 0.08065034 
... Procrustes: rmse 2.805412e-05  max resid 5.594612e-05 
... Similar to previous best
Run 3 stress 0.1378308 
Run 4 stress 0.08065035 
... Procrustes: rmse 7.802029e-05  max resid 0.0001548812 
... Similar to previous best
Run 5 stress 0.08065034 
... Procrustes: rmse 3.731121e-05  max resid 6.566436e-05 
... Similar to previous best
Run 6 stress 0.1378308 
Run 7 stress 0.08554815 
Run 8 stress 0.08065037 
... Procrustes: rmse 0.0001137107  max resid 0.0002575052 
... Similar to previous best
Run 9 stress 0.08554815 
Run 10 stress 0.1378309 
Run 11 stress 0.08065034 
... New best solution
... Procrustes: rmse 1.388738e-05  max resid 1.990997e-05 
... Similar to previous best
Run 12 stress 0.1378309 
Run 13 stress 0.1905268 
Run 14 stress 0.1743257 
Run 15 stress 0.08065036 
... Procrustes: rmse 7.179091e-05  max resid 0.0001453003 
... Similar to previous best
Run 16 stress 0.08065034 
... Procrustes: rmse 4.443671e-05  max resid 7.380012e-05 
... Similar to previous best
Run 17 stress 0.08065034 
... Procrustes: rmse 9.112915e-06  max resid 1.854915e-05 
... Similar to previous best
Run 18 stress 0.08065035 
... Procrustes: rmse 4.869656e-05  max resid 9.470507e-05 
... Similar to previous best
Run 19 stress 0.1460913 
Run 20 stress 0.1378308 
*** Best solution repeated 5 times

nmds.2021.2

Call:
metaMDS(comm = matrix.bray.2021.2, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021.2 
Distance: bray 

Dimensions: 2 
Stress:     0.08065034 
Stress type 1, weak ties
Best solution was repeated 5 times in 20 tries
The best solution was from try 11 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021.2’ 

stressplot(nmds.2021.2)

plot(nmds.2021.2, type="t")

dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))

nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))

nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")

nmds 2021

Stress: fair

0.1230831

metaMDS()

nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)

Run 0 stress 0.1232098 
Run 1 stress 0.1339879 
Run 2 stress 0.1230831 
... New best solution
... Procrustes: rmse 0.04463789  max resid 0.1321469 
Run 3 stress 0.1230831 
... Procrustes: rmse 1.350077e-05  max resid 4.219102e-05 
... Similar to previous best
Run 4 stress 0.133988 
Run 5 stress 0.1507047 
Run 6 stress 0.1573377 
Run 7 stress 0.1230831 
... Procrustes: rmse 2.810056e-06  max resid 8.291934e-06 
... Similar to previous best
Run 8 stress 0.1339879 
Run 9 stress 0.1230831 
... Procrustes: rmse 8.646086e-06  max resid 2.745517e-05 
... Similar to previous best
Run 10 stress 0.132972 
Run 11 stress 0.1232098 
... Procrustes: rmse 0.04463824  max resid 0.1315186 
Run 12 stress 0.1230831 
... New best solution
... Procrustes: rmse 2.871318e-06  max resid 7.467703e-06 
... Similar to previous best
Run 13 stress 0.1421227 
Run 14 stress 0.1232098 
... Procrustes: rmse 0.04463976  max resid 0.1315202 
Run 15 stress 0.1339879 
Run 16 stress 0.1515614 
Run 17 stress 0.1230832 
... Procrustes: rmse 2.468942e-05  max resid 7.756011e-05 
... Similar to previous best
Run 18 stress 0.1573377 
Run 19 stress 0.1230831 
... Procrustes: rmse 7.42768e-06  max resid 2.388952e-05 
... Similar to previous best
Run 20 stress 0.132972 
*** Best solution repeated 3 times

nmds.2021

Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021 
Distance: bray 

Dimensions: 2 
Stress:     0.1230831 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 12 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021’ 

stressplot(nmds.2021)

plot(nmds.2021, type="t")

nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
  )
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))

nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Precipitation Treatment")

2021 — blocks

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block <- matrix0.block %>%
  filter(yr == 2021)

row.2021.block <- paste(matrix.2021.block$nut_trt, matrix.2021.block$ppt_trt, matrix.2021.block$yr, matrix.2021.block$block, sep = "_")

matrix.2021.block <- matrix.2021.block[, c(5:79)]

rownames(matrix.2021.block) <- row.2021.block

Warning: Setting row names on a tibble is deprecated.

matrix.bray.block.2021 <- decostand(matrix.2021.block, "total")

nmds 2021

Stress: fair

0.1702944

metaMDS()

nmds.2021.block <- metaMDS(matrix.bray.block.2021, k=2, trymax = 25)

Run 0 stress 0.1702944 
Run 1 stress 0.2323726 
Run 2 stress 0.2333131 
Run 3 stress 0.1702944 
... New best solution
... Procrustes: rmse 8.12052e-06  max resid 2.925004e-05 
... Similar to previous best
Run 4 stress 0.1704654 
... Procrustes: rmse 0.006480506  max resid 0.02760662 
Run 5 stress 0.2085283 
Run 6 stress 0.2044772 
Run 7 stress 0.2146073 
Run 8 stress 0.1942847 
Run 9 stress 0.1873962 
Run 10 stress 0.2098049 
Run 11 stress 0.1917675 
Run 12 stress 0.1833145 
Run 13 stress 0.2044772 
Run 14 stress 0.226506 
Run 15 stress 0.2142311 
Run 16 stress 0.1702944 
... Procrustes: rmse 5.781581e-06  max resid 1.825294e-05 
... Similar to previous best
Run 17 stress 0.1798315 
Run 18 stress 0.1704654 
... Procrustes: rmse 0.006481124  max resid 0.02760841 
Run 19 stress 0.2355275 
Run 20 stress 0.2089454 
*** Best solution repeated 2 times

nmds.2021.block

Call:
metaMDS(comm = matrix.bray.block.2021, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021 
Distance: bray 

Dimensions: 2 
Stress:     0.1702944 
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 3 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021’ 

stressplot(nmds.2021.block)

plot(nmds.2021.block, type="t")

nmds1.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 4)
  )
nmds2.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(2), display=c("sites")))

nmds_2021_dat.block <- cbind(nmds1.2021.block, nmds2.2021.block) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Precipitation Treatment")

2021 — blocks (LOW)

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
  filter(yr == 2021,
         block == 3 | block == 4)

row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")

matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]

rownames(matrix.2021.block.low) <- row.2021.block.low

Warning: Setting row names on a tibble is deprecated.

matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")

nmds 2021

Stress: fair

0.1282296

metaMDS()

nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)

Run 0 stress 0.1282296 
Run 1 stress 0.1282296 
... New best solution
... Procrustes: rmse 6.900845e-05  max resid 0.0001846812 
... Similar to previous best
Run 2 stress 0.1282296 
... New best solution
... Procrustes: rmse 6.471065e-05  max resid 0.0002196127 
... Similar to previous best
Run 3 stress 0.1282296 
... Procrustes: rmse 8.428118e-05  max resid 0.0003017723 
... Similar to previous best
Run 4 stress 0.1985506 
Run 5 stress 0.1282296 
... Procrustes: rmse 1.033532e-05  max resid 2.797048e-05 
... Similar to previous best
Run 6 stress 0.1785638 
Run 7 stress 0.1282296 
... Procrustes: rmse 8.012815e-05  max resid 0.0002839093 
... Similar to previous best
Run 8 stress 0.1282296 
... Procrustes: rmse 5.434575e-05  max resid 0.0001894526 
... Similar to previous best
Run 9 stress 0.1282296 
... New best solution
... Procrustes: rmse 1.600545e-05  max resid 5.619322e-05 
... Similar to previous best
Run 10 stress 0.1282296 
... Procrustes: rmse 1.523127e-05  max resid 5.391544e-05 
... Similar to previous best
Run 11 stress 0.1428149 
Run 12 stress 0.1627541 
Run 13 stress 0.2101121 
Run 14 stress 0.1282296 
... Procrustes: rmse 2.043501e-05  max resid 6.96312e-05 
... Similar to previous best
Run 15 stress 0.2053314 
Run 16 stress 0.2106634 
Run 17 stress 0.1282296 
... Procrustes: rmse 8.924757e-05  max resid 0.0003091915 
... Similar to previous best
Run 18 stress 0.1428149 
Run 19 stress 0.1428149 
Run 20 stress 0.1428149 
*** Best solution repeated 4 times

nmds.2021.block.low

Call:
metaMDS(comm = matrix.bray.block.2021.low, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.low 
Distance: bray 

Dimensions: 2 
Stress:     0.1282296 
Stress type 1, weak ties
Best solution was repeated 4 times in 20 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.low’ 

stressplot(nmds.2021.block.low)

plot(nmds.2021.block.low, type="t")

nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
  )
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))

nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")

2021 — blocks (HIGH)

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.high <- matrix0.block %>%
  filter(yr == 2021,
         block == 1 | block == 2)

row.2021.block.high <- paste(matrix.2021.block.high$nut_trt, matrix.2021.block.high$ppt_trt, matrix.2021.block.high$yr, matrix.2021.block.high$block, sep = "_")

matrix.2021.block.high <- matrix.2021.block.high[, -c(0:4)]

rownames(matrix.2021.block.high) <- row.2021.block.high

Warning: Setting row names on a tibble is deprecated.

matrix.bray.block.2021.high <- decostand(matrix.2021.block.high, "total")

nmds 2021

Stress: fair

0.1412997

metaMDS()

nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)

Run 0 stress 0.1412997 
Run 1 stress 0.1430823 
Run 2 stress 0.1430824 
Run 3 stress 0.1856239 
Run 4 stress 0.1581575 
Run 5 stress 0.1412999 
... Procrustes: rmse 0.0003384719  max resid 0.001097121 
... Similar to previous best
Run 6 stress 0.1430825 
Run 7 stress 0.1581574 
Run 8 stress 0.1856239 
Run 9 stress 0.1412997 
... Procrustes: rmse 3.643604e-05  max resid 0.0001217465 
... Similar to previous best
Run 10 stress 0.1581574 
Run 11 stress 0.1581574 
Run 12 stress 0.1430824 
Run 13 stress 0.1430824 
Run 14 stress 0.1581574 
Run 15 stress 0.2412356 
Run 16 stress 0.1430827 
Run 17 stress 0.1430824 
Run 18 stress 0.1412998 
... Procrustes: rmse 0.0002267585  max resid 0.0007337562 
... Similar to previous best
Run 19 stress 0.1581575 
Run 20 stress 0.1581574 
*** Best solution repeated 3 times

nmds.2021.block.high

Call:
metaMDS(comm = matrix.bray.block.2021.high, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.high 
Distance: bray 

Dimensions: 2 
Stress:     0.1412997 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.high’ 

stressplot(nmds.2021.block.high)

plot(nmds.2021.block.high, type="t")

nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
  )
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))

nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")

Make species x trait matrix

Greenhouse traits

• height, leaf dry matter content, specific leaf area, root mass fraction, relative growth rate, coarse root diameter, root density, specific root length of fine roots, specific root length of coarse roots

Field traits

• height, leaf dry matter content, specific leaf area, seed mass, leaf C, leaf N

“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”

“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”

Workflow plan:

• overlay species trait data as vectors on the NMDS

  • only greenhouse traits available (only sla data from field)

    • use PC scores as composite trait axes OR

    • average the traits for each species

Read in trait data

# # read in greenhouse trait data
# traits.gh <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=ya1ck6p7&dl=1")
# 
# # read in field trait data
# # download.file("https://www.dropbox.com/scl/fi/trqrnz54dkjoa9do4tle5/leaf-area-data_NG.xlsx?rlkey=s0kic2a5yqknikc4icek8y7rn&st=dl4x634h&dl=1",destfile = "field_sla.xlsx")
# traits.sla.field <- 
#   read_excel("field_sla.xlsx",col_names = TRUE,col_types = "text")
# 
# rda(traits.gh, scale = TRUE)

Investigate trait data, looking in Greenhouse trait folder

Trait data is identical. Will clean data in

# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
# 
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
# 
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)

traits.cleaned <- read.csv("data/GreenhouseTraits_corrected_20240212.csv") %>%
  mutate(
    Fresh.leaf.mass..g. = as.numeric(ifelse(
      Fresh.leaf.mass..g. == "Skipped accidentally", NA, Fresh.leaf.mass..g.
    )),
    Dry.leaf.mass..g. = as.numeric(ifelse(
      Dry.leaf.mass..g. == "No sample", NA, Dry.leaf.mass..g.
    )),
    Root.dry.biomass..g. = as.numeric(ifelse(
      Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.
    )),
    Root.volume..cm3. = as.numeric(ifelse(
      Root.volume..cm3. == "not scanned", NA, Root.volume..cm3.
    ))
  ) 

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Root.dry.biomass..g. = as.numeric(ifelse(Root.dry.biomass..g. == "roots lost", NA,
  Root.dry.biomass..g.))`.
Caused by warning:
! NAs introduced by coercion

Summarize trait data (GH) for each species

means

# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
  mutate(
    ID = ifelse(
      ID == "Agoseris", "AGSP", ifelse(
        ID == "AVEBAR", "AVBA", ifelse(
          ID == "BRDR", "BRDI", ID
        )
      )
    )
  )

# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
  group_by(ID) %>%
  summarise(
    across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
  )))
  

# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)

intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)

 [1] "ASSP"  "AST1"  "AST2"  "AST3"  "CIQU"  "CRTI"  "CYEC"  "DACA"  "DIMU"  "GAGR"  "GAPA"  "GAPH"  "HYSP" 
[14] "JUBU"  "LYHY"  "MAD1"  "MAD2"  "NAPU"  "PAVI"  "SABI1" "SABI2" "SIGA"  "THGR"  "TRSP"  "UNBU"  "UNF1" 
[27] "UNF3"  "UNF4"  "UNF8"  "UNFR"  "UNGR"  "XAST" 

setdiff(traits.means$ID, colnames(matrix.bray))

 [1] "ACMI"      "ANARp"     "AVFA"      "BRCA"      "BRHOp"     "BRNIf"     "BRNIp"     "CLPUf"     "CLPUp"    
[10] "CYDA"      "Crassula"  "ESCA"      "FEMI"      "FEMY"      "GITRf"     "GITRp"     "HYGL"      "LACA"     
[19] "LENIf"     "LENIp"     "LOPU"      "LUBI"      "Leontodon" "MAELf"     "MAELp"     "MICA"      "NEMA"     
[28] "PLERf"     "PLERp"     "RUPU"      "Sanicula"  "TRCI"      "TRHIpi"    "TRWIf"     "TRWIp"     "TRWIpi"   
[37] "VIVA"     

# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
  dplyr::select(code4, species) %>%
  summarise(
    code4 = unique(code4),
    species = unique(species)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.

unique.trait.species <- traits.cleaned %>%
  dplyr::select(Taxon, ID) %>%
  distinct(ID, .keep_all = TRUE)
  

• change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)

• What’s the difference between:

  • ‘ANAR’ and ‘ANARp’ in trait data?

  • ‘BRHO’ and ‘BRHOp’ in trait data?

  • ‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)

  • ‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)

  • ‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)

  • ‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)

  • ‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)

  • ‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)

  • ‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)

  • ‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)

  • ‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)

  • ‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)

  • ‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)

• change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)

• change BRDR –> BRDI

Overlay Trait Vectors onto NMDS

ALL YEARS – LOW AND HIGH GRAZING

Prepare NMDS and Traits

Stress:

0.2011127

# subset nmds to match species names
nmds.comp.trait <- metaMDS(matrix.bray %>%
                             dplyr::select(intersect.names), 
                           k=2, trymax = 25)

Run 0 stress 0.2060005 
Run 1 stress 0.2007525 
... New best solution
... Procrustes: rmse 0.06486829  max resid 0.1521338 
Run 2 stress 0.2133228 
Run 3 stress 0.2079153 
Run 4 stress 0.2051687 
Run 5 stress 0.2107849 
Run 6 stress 0.2078775 
Run 7 stress 0.2027188 
Run 8 stress 0.2013372 
Run 9 stress 0.2050098 
Run 10 stress 0.2017137 
Run 11 stress 0.2042866 
Run 12 stress 0.2051682 
Run 13 stress 0.2185014 
Run 14 stress 0.2078775 
Run 15 stress 0.2156872 
Run 16 stress 0.2221242 
Run 17 stress 0.2079153 
Run 18 stress 0.2023083 
Run 19 stress 0.2181406 
Run 20 stress 0.2010909 
... Procrustes: rmse 0.08100757  max resid 0.2166622 
Run 21 stress 0.2086475 
Run 22 stress 0.2046958 
Run 23 stress 0.2061378 
Run 24 stress 0.2020872 
Run 25 stress 0.2080994 
*** Best solution was not repeated -- monoMDS stopping criteria:
    24: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

nmds.comp.trait # pretty not great stress

Call:
metaMDS(comm = matrix.bray %>% dplyr::select(intersect.names),      k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray %>% dplyr::select(intersect.names) 
Distance: bray 

Dimensions: 2 
Stress:     0.2007525 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray %>% dplyr::select(intersect.names)’ 

stressplot(nmds.comp.trait)

plot(nmds.comp.trait, type="t")

nmds1 <- as.data.frame(scores(nmds.comp.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.comp.trait, choices=c(2), display=c("sites")))

nmds_dat_trait <- cbind(nmds1, nmds2) %>%
  as.data.frame() %>%  # Ensure it's a data frame
  rownames_to_column(var = "site_id") %>%  # Move row names into a new column
  separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")


# subset matrix and trait data to match names
traits.means.comp <- traits.means %>%
  filter(ID %in% intersect.names) %>%
  rename(
    Height = Height..cm._mean,
    FreshLeafMass = Fresh.leaf.mass..g._mean,
    DryLeafMass = Dry.leaf.mass..g._mean,
    LDMC = LDMC_mean, # Leaf dry matter content (LDMC, the ratio of leaf dry mass to fresh mass)
    LeafArea = Leaf.Area..cm2._mean,
    SLA = SLA..cm2.g._mean,
    ShootDryBiomass = Shoot.dry.biomass..g._mean,
    RootDryBiomass = Root.dry.biomass..g._mean,
    TotalBiomass = Total.biomass..g._mean,
    RMF = RMF_mean, # Root mass fraction (RMF) is a plant trait that measures the proportion of a plant's dry mass that is in its roots
    RootVol = Root.volume..cm3._mean,
    RootDensity = Root.density..g.cm3._mean,
    CoarseRootDiameter = Coarse.root.diameter..mm._mean,
    RootLength = Length..mm._mean,
    FineRootLength = Fine.root.length..mm._mean,
    CoarseRootLength = Coarse.root.length..mm._mean,
    CoarseRootSpecLength = Coarse.root.specific.length..cm.g._mean,
    FineRootSpecLength = Fine.root.specific.length..cm.g._mean,
    PropFineRoots = Proportion.fine.roots_mean
    
  )
# traits.means.comp2 <- traits.means.comp %>%
# dplyr::select(-ID)
rownames(traits.means.comp2) <- traits.means.comp$ID

Warning: Setting row names on a tibble is deprecated.

# trait vectors onto the NMDS
# envfit(nmds.comp.trait, traits.means.comp, permutations = 999)
# example from stack overflow: envfit(ord ~ var1 + var2 + var3, dune.spec, display="sp"), 'ord' is the nmds model, 'var1' is the trait1 i think?

colnames(traits.means.comp)

 [1] "ID"                   "Height"               "FreshLeafMass"        "DryLeafMass"         
 [5] "LDMC"                 "LeafArea"             "SLA"                  "ShootDryBiomass"     
 [9] "RootDryBiomass"       "TotalBiomass"         "RMF"                  "RootVol"             
[13] "RootDensity"          "CoarseRootDiameter"   "RootLength"           "FineRootLength"      
[17] "CoarseRootLength"     "CoarseRootSpecLength" "FineRootSpecLength"   "PropFineRoots"       

vectors.nmds.comp.trait <- envfit(nmds.comp.trait ~ Height + 
                                    LeafArea +
                                    RMF, traits.means.comp, display = "sp") # ONLY Height, Leaf area, and RMF significant

# is it different when you only do certain traits (above- vs. below-ground)? or can you just throw everything in there and it will tell you what is significant and you can subset from there?

# Looks like different subsets of traits put in are the same, so I shall proceed with the mega models and verify with Lauren

trait_vectors <- as.data.frame(vectors.nmds.comp.trait$vectors$arrows)

NMDS model: only Height, Leaf area, and RMF () found significant.

Plot

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  # stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
    values = c(
      "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
      "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
      "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
      "wet" = rgb(0, 0, 255, maxColorValue = 255),
      "dry" = rgb(183, 65, 14, maxColorValue = 255)
    )
  ) +
  geom_segment(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
    inherit.aes = FALSE,  # 🚀 Prevents it from looking for 'nut_trt'
    arrow = arrow(length = unit(0.2, "cm")), 
    color = "red", 
    size = 1
  ) +
  geom_text(
    data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows), 
    aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)), 
    inherit.aes = FALSE,  # Avoids unwanted inheritance
    hjust = -0.2, 
    vjust = -0.2, 
    color = "red",
    size = 5
  )

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)